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1 Introduction 

The flow shop scheduling problem consists in the assignment of a set of jobs J = {Ji, . . . , J n }, 
each of which consists of a set of operations Jj = {Oji, . . . ,Oj Qj } onto a set of machines M = 
{Mi, . . . , M m } [5, 18]. Each operation Ojk is processed by at most one machine at a time, involving 
a non-negative processing time pjk- The result of the problem resolution is a schedule x, defining for 
each operation Ojk a starting time Sjk on the corresponding machine. Several side constraints are 
present which have to be respected by any solution x belonging to the set of feasible schedules X. 
Precedence constraints Ojk t> Ojfc+iVj = 1, . . . , n, k = 1, . . . , Oj — 1 between the operations of a job 
Jj assure that processing of 0jk+i only commences after completion of Ojk, thus Sjk+i ^ Sjk +Pjk- 
In flow shop scheduling, the machine sequence in which the operations are processed by the machines 
is identical for all jobs, and for the specific case of the permutation flow shop scheduling the job 
sequence must also be the same on all machines. 

The assignment of starting times to the operations has to be done with respect to one or several 
optimality criteria. Most optimality criteria are functions of the completion times Cj of the jobs 

J 3i = s joj + PjOj 

The most prominent optimality criteria are the maximum completion time (makespan) C max = 
rn.ax.Cj and the sum of the completion times C sum = Y^j=\^j- Others express violations of due 
dates dj of jobs Jj. A due date dj defines a latest point of time until production of a job Jj should 
be completed as the finished product has to be delivered to the customer on or up to this date. A 
possible optimality criteria based on tardiness of jobs is e.g. the total tardiness T sum = X]y=i^}> 
where Tj = max(Cj — dj,0). 

It is known, that for regular optimality criteria at least one active schedule x does exist which 
is also optimal. The representation of an active schedule is possible using a permutation of jobs 
7r = (-7Ti, . . . , 7r n ), where each Tij stores a job Jk at position j. The permutation is then decoded into 
an active schedule by assuming the starting times of all operations as early as possible with respect 
to the precedence constraints and the given sequence in ir. As a consequence, the search is, instead 
of searching all possible schedules, restricted to the much smaller set of active schedules only. 
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Multi-objective approaches to scheduling consider a vector G(x) = (gi(x), . . . ,gic(x)) of opti- 
mally criteria at once [21]. As the relevant optimality criteria are often of conflicting nature, not 
a single solution x € X exists optimizing all components of G(x). Optimality in multi-objective 
optimization problems is therefore understood in the sense of Pareto-optimality, and the resolution 
of multi-objective optimization problems lies in the identification of all elements belonging to the 
Pareto set P, containing all alternatives x which are not dominated by any other alternative x' £ X. 

Several approaches of metaheuristics have been formulated and tested in order to solve the per- 
mutation flow shop scheduling problem under multiple, in most cases two, objectives. Common to 
all is the representation of solutions using permutations ir of jobs, as in previous investigations only 
regular functions are considered. 

First results have been obtained using Evolutionary Algorithms, which in general play a dominant 
role in the resolution of multi-objective optimization problems when using metaheuristics. This is 
mainly due to the fact that these methods incorporate the idea of a set of solutions, a so called 
population, as a general ingredient. Flow shop scheduling problems minimizing the maximum com- 
pletion time and the average flow time have been solved by Nagar, Heragu and Haddock [17]. In 
their work, they however combine the two objectives into a weighted sum. Problems minimizing the 
maximum completion time and the total tardiness are solved by Murata, Ishibuchi and Tanaka 
[16], again under the combination of both objectives into a weighted sum. Later work on the same 
problem class by Basseur, Seynhaeve and Talbi [4] avoids the weighted sum approach, using 
dominance relations among the solutions only. 

Most recent work is presented by Loukil, Teghem and Tuyttens [14]. Contrary to approaches 
from Evolutionary Computations, the authors apply the Multi Objective Simulated Annealing ap- 
proach MOSA of Ulungu, Teghem, Fortemps and Tuyttens [22] to a variety of bi-criterion 
scheduling problems. 

Flow shop scheduling problems with three objectives are studied by Ishibuchi and Murata [11], 
and Ishibuchi, Yoshida and Murata [12]. The authors minimize the maximum completion time, 
the total completion time, and the maximum tardiness at once. A similar problem minimizing the 
maximum completion time, the average flow time, and the average tardiness is then tackled by 
Bagchi [1, 2]. 



2 Pareto Iterated Local Search 



The Pareto Iterated Local Search (PILS) metaheuristic is a concept for the solution of multi- 
objective optimization problems. It combines the two main driving forces of local search, inten- 
sification and diversification, into a single algorithm, and extends the work presented in [8]. The 
motivation behind this concept can be seen in the increasing demand for simple, yet effective heuris- 
tics for the resolution of complex multi-objective optimization problems. Two developments in local 
search demonstrate the effectiveness of some intelligent ideas that make use of certain structures 
within the search space topology of problems: 

1. Iterated Local Search, introducing the idea of perturbing solutions to overcome local optimality 
and continue search in interesting areas of the search space [15]. After the pioneering work of 
Boese [6], who investigated properties of the search space of the traveling salesman problem, 
this concept has been used with increasing success on problems where solutions of high quality 
can be found relatively concentrated in alternative space. 
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2. Second, Variable Neighborhood Search [10], combining multiple neighborhood operators into 
a single algorithm in order to avoid local optimality in the first place. 



In the proposed concept, both paradigms are combined and extended within a search framework 
handling not only a single but a set of alternatives at once. 

The main principle of the algorithm is sketched in Figure 1. Starting from an initial solution 
xi, an improving, intensifying search is performed until a set of locally optimal alternatives is 
identified, stored in a set p a PP rox representing the approximation of the true Pareto set P. No 
further improvements are possible from this point. In this initial step, a set of neighborhoods ensures 
that all identified alternatives are locally optimal not only to a single but to a set of neighborhoods. 
This principle, known from Variable Neighborhood Search, promises to lead to better results as it 
is known that all global optima are also locally optimal with respect to all possible neighborhoods 
while this is not necessarily the case for local optima. 



G(x,)o 




Figure 1: Illustration of the Pareto Iterated Local Search metaheuristic. The archive of the currently 
best solutions is updated during the search. Here, G{x^) dominates G{x2) which is going to be 
deleted from p a PP rox . 



After the identification of a locally optimal set, a diversification step is performed on a solution 
X2 using a perturbation operator, continuing search from the perturbed solution X3. The pertur- 
bation operator has to be significantly different from the neighborhoods used in intensification, as 
otherwise the following search would return to the previous solution. On the other hand however, 
the perturbation should not entirely destroy the characteristics of the alternative. Doing that would 
result in a random restart of the search without keeping promising attributes of solutions. 

The PILS metaheuristic may be formalized as given in Algorithm 1. The intensification of the 
algorithm, illustrated in steps (1) and (3) of Figure 1 is between the lines 6 to 21, the description 
of the diversification, given in step (2) of Figure 1 is between the lines 22 to 26. 

It can be seen, that the algorithm computes a set of neighborhoods for each alternative. The 
sequence in which the neighborhoods are computed is arranged in a random fashion, described 
in line 13 of Algorithm 1. This introduces an additional element of diversity to the concept, as 
otherwise the search might be biased by a certain sequence of neighborhoods. 
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Algorithm 1 Pareto Iterated Local Search 
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Set x = x' 
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until x locally optimal with respect to Nx, . . . , N^ 


, therefore i > k 
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Set neighborhoods of x as 'investigated' 
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Set i = 1 




on 
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if 3x' £ p a PP rox neighborhoods not investigated yet then 


21 


Set x = x' 
else 




22 




23 


Randomly select some x' 6 papprox 




24 


Compute x" = N perturb (x'\ 




25 


Set x = x" 




26 


end if 




27 


until termination criterion is met 







3 Experiments on mult i- objective flow shop scheduling problems 



3.1 Algorithm configuration and experimental setup 

In the following, the Pareto Iterated Local Search is applied to a set of benchmark instances of the 
multi-objective permutation flow shop scheduling problem. The first instances have been provided 
by Basseur, Seynhaeve and Talbi [4], who defined due dates for the well-known instances of 
Taillard [20]. The instances range from n = 20 jobs that have to be processed on m = 5 
machines to n = 100, m = 20. All of them are solved under the simultaneous consideration of the 
minimization of the maximum completion time C max and the total tardiness T sum and are referred 
to as 'Ta n x m'. 

We also solved the benchmark instance 'Ba 49 x 15' by Bagchi [1], consisting of n = 49 jobs 
on m = 15 machines. The three objective functions of the data set are the minimization of the 
maximum completion time C max , the minimization of the average completion time ^C sum , and the 
minimization of the average tardiness z:T surri . 

Three operators are used in the definition of the neighborhoods Nx, . . . ,Nfc, described in the 
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work of Reeves [19]. First, an exchange neighborhood, exchanging the position of two jobs in it, 
second, a forward shift neighborhood, taking a job from position i and reinserting it at position 
j with j < i, and finally a backward shift neighborhood, shifting a job from position i to j with 
j < i. All operators are problem independent operators for permutation-based representations, each 
computing n ^ n 2 ^ neighboring solutions. 

After a first approximation p a PP rox of the Pareto set is obtained, one element x' p a PP rox is 
selected by random and perturbed into another solution x". We use a special neighborhood that on 
one hand leaves most of the characteristics of the perturbed alternatives intact, while still changes 
the positions of some jobs. Also, several consecutive applications of the neighborhoods Ni, . . . , 
would be needed to transform x" back into x' . This is important, as otherwise the algorithm might 
return to the initially perturbed alternative x', leading to a cycle in the search path. 
The perturbation neighborhood N per turb can be described as follows. First, a subset of ir is randomly 
selected, comprising four consecutive jobs at positions j, j+l,j+2, j+3. Then a neighboring solution 
x" is generated by moving the job at position j to j + 3, the one at position j + 1 to j + 2, the one 
at position j + 2 to j, and the job at position j + 3 to j + 1, leaving the jobs at the positions < j 
and > j + 3 untouched. In brief, this leads to a combination of several exchange and shift moves, 
executed at once. 

In order to analyze the quality of the approximations, we compare the results obtained by 
PILS to the approximations of a multi-objective multi-operator search algorithm MOS, described 



in 


Algorithm 2. 




Algorithm 2 Multi-objective multi-operator search framework 
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Generate initial solution x, set p a PP rox = {x} 
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repeat 




3 


Randomly select some x € papprox neighborhoods 


not investigated yet 


4 


Randomly select some neighborhood Nj from Ni, . 




r> 


Generate Nj(x) 
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Update p a PP rox with N^x) 




7 


it x £ p a PP r ° x then 




<S 


Set neighborhoods of x as 'investigated' 




9 


end if 




10 


until ~fix € p a PP rox | neighborhoods not investigated yet 


11 


Return P a PP r ° x 





The MOS Algorithm, taken from [8], is based on the concept of Variable Neighborhood Search, 
extending the general idea of several neighborhood operators by adding an archive p a PP rox towards 
the optimization of multi-objective problems. For a fair comparison, the same neighborhood opera- 
tors are used as in the PILS algorithm. After the termination criterion is met in step 10, we restart 
search while keeping the approximation p a w rox for the final analysis of the quality of the obtained 
solutions. 

3.2 Results 

The benchmark instances of Basseur and the one of Bagchi have been solved using the PILS 
algorithm. In each of the 100 test runs, the approximation quality of the obtained results has been 
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analyzed using the D\ and Z?2 metrics of Czyzak and Jaszkiewicz [7]. The two metrics have been 
chosen for the analysis as they provide an interesting interpretation from an economical point of 
view. Based on a so called 'achievement scalarizing function', they compute the average {D\) and 
the maximum (-D2) regret a decision maker would have to face when trying to select a certain most 
preferred alternative x* G P, approximated by the results in p a PP rox m 

While for the smaller instances the optimal solutions are known, the analysis for the larger 
instances has to rely on the best known results published in the literature. Experiments have been 
carried out on a Intel Pentium IV processor, running at 1.8 GHz. Table 1 gives an overview about 
the number of evaluations executed for each instance. Clearly, considerable more alternatives have 
to be evaluated with increasing size of the problem instances to allow a convergence of the algorithm. 
Also, the running times, given in Table 1, too, increase with increasing size of the problem instances. 
No significant difference in running behavior can be found when comparing the two metaheuristics 
PILS and MOS. Apart from some minor differences around the perturbation neighborhood N per t ur &, 
the approaches are identical with respect to the impact on the resulting running times as they use 
the same neighborhood operators. 

Table 1: Number of evaluations and running times for the investigated instances. 



Instance n x m No of evaluations Eval. time Neighbor. 

comp. time 



Ta 20 x 5 (#1) 


1,000,000 


42.7 


0.5 


Ta 20 x 5 (#2) 


1,000,000 


42.5 


0.6 


Ta 20 x 10 (#1) 


1,000,000 


78.0 


0.6 


Ta 20 x 10 (#2) 


1,000,000 


83.3 


0.6 


Ta 20 x 20 


1,000,000 


143.5 


0.6 


Ta 50 x 5 


5,000,000 


195.1 


1.3 


Ta 50 x 10 


5,000,000 


386.1 


1.3 


Ta 50 x 20 


5,000,000 


754.3 


1.3 


Ta 100 x 10 


10,000,000 


1459.7 


2.5 


Ta 100 x 20 


10,000,000 


2885.3 


2.5 


Ba 49 x 15 


5,000,000 


566.3 


1.2 



Times are given in milliseconds. 



An implementation of the algorithm has been made available within an integrated software for 
the resolution of multi-objective scheduling problems using metaheuristics. The system is equipped 
with an extensive user interface that allows an interaction with a decision maker and is able to visu- 
alize the obtained results in alternative and outcome space. The system also allows the comparison 
of results obtained by different metaheuristics. 

The average values obtained by the investigated metaheuristics are given in Table 2. It can be 
seen, that PILS leads for all investigated problem instances to better results for both the Di and 
the D2 metric. This general result is consistently independent from the actual problem instance and 
significant at a level of significance of 0.01. For a single instance, the 'Ta 20 x 5 (#1)', PILS was 
able to identify all optimal solutions in all test runs, leading to average values of D\ = D2 = 0.0000. 
Apparently, this instance is comparably easy to solve. 
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Table 2: Average results of D\ and D 2 



Di D 2 



Instance n x m 


PILS 


MOS 


PILS 


MOS 


Ta 20 x 5 (#1) 


0.0000 


0.0323 


0.0000 


0.1258 


Ta 20 x 5 (#2) 


0.1106 


0.1372 


0.3667 


0.4249 


Ta 20 x 10 (#1) 


0.0016 


0.0199 


0.0146 


0.0598 


Ta 20 x 10 (#2) 


0.0011 


0.0254 


0.0145 


0.1078 


Ta 20 x 20 


0.0088 


0.0286 


0.0400 


0.1215 


Ta 50 x 5 


0.0069 


0.0622 


0.0204 


0.1119 


Ta 50 x 10 


0.0227 


0.3171 


0.0897 


0.4658 


Ta 50 x 20 


0.0191 


0.3966 


0.0616 


0.5609 


Ta 100 x 10 


0.0698 


0.3190 


0.1546 


0.4183 


Ta 100 x 20 


0.0013 


0.2349 


0.0255 


0.3814 


Ba 49 x 15 


0.0202 


0.2440 


0.0701 


0.3414 



While it was possible to show in [8] that the MOS algorithm is competitive to different Evolu- 
tionary Algorithms, iterating search in qualitatively good areas of the search space by PILS improves 
the results even further. Recalling, that with increasing problem size an increasing amount of time 
is needed to evaluate the alternatives, iterating in promising regions becomes even more interesting 
as opposed to restarting the search. 

A deeper analysis has been performed to monitor the resolution behavior of the local search 
algorithms and to get a better understanding of how the algorithm converges towards the Pareto 
front. Figure 2 (a) plots with symbol x the results obtained by random sampling 50,000 alter- 
natives for the problem instance 'Ta 100 x 10', and compares the points obtained during the first 
intensification procedure of PILS until a locally optimal set is identified. The alternatives computed 
starting from a random initial solution towards the Pareto front are plotted as +, the Pareto front 
as 0. It can be seen, that in comparison to the initial solution even a simple local search approach 
converges in rather close proximity to the Pareto front. With increasing number of computations 
however, the steps towards the optimal solutions get increasingly smaller, as it can be seen when 
monitoring the distances between the + symbols. After convergence towards a locally optimal set, 
overcoming local optimality is then provided by means of the perturbation neighborhood N per t ur fe. 

An interesting picture is obtained when analyzing the distribution of the randomly sampled 
50,000 alternatives for instance Ta 100 x 10'. In Figure 2 (b), the number of alternatives with 
a certain combination of objective function values are plotted and compared to the Pareto front, 
given in the left corner. It turns out that many alternatives are concentrated around some value 
combination in the area of approximately C max = 6900, T sum = 111500, relatively far away from 
the Pareto front. 

When analyzing the convergence of local search heuristics toward the globally Pareto front 
as well as towards locally optimal alternatives, the question arises how many local search steps 
are necessary until a locally optimal alternative is identified. From a different point of view, this 
problem is discussed in the context of computational complexity of local search [13]. It might 
be worth investigating this behavior in quantitative terms. Table 3 gives the average number of 
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Figure 2: (a) Randomly generated solutions (x), intensification of search(+), and Pareto front (©). 
(b) Distribution of randomly generated solutions (+) compared to the Pareto front (©) 

evaluations that have been necessary to reach a locally optimal alternative from some randomly 
generated initial solution. The analysis reveals that the computational effort grows exponentially 
with the number of jobs n. 



Table 3: Average number of evaluations until a locally optimal alternative is reached 



Instance n x m 


No of jobs 


No of evaluations 


Instance n x m 


No of jobs 


No of eval. 


Ta 20 x 5 (#1) 


20 


3,614 


Ta 50 x 5 


50 


53,645 


Ta 20 x 5 (#2) 


20 


3,292 


Ta 50 x 10 


50 


55,647 


Ta 20 x 10 (#1) 


20 


2,548 


Ta 50 x 20 


50 


38,391 


Ta 20 x 10 (#2) 


20 


2,467 


Ta 100 x 10 


100 


793,968 


Ta 20 x 20 


20 


2,657 


Ta 100 x 20 


100 


479,420 








Ba 49 x 15 


49 


28,908 



4 Conclusions 

In the past years, considerable progress has been made in solving complex multi-objective optimiza- 
tion problems. Effective metaheuristics have been developed, providing the possibility of computing 
approximations to problems with numerous objectives and complex side constraints. While many 
approaches are of increasingly effectiveness, complex parameter settings are however required to 
tune the solution approach to the given problem at hand. 

The algorithm presented in this paper proposed a metaheuristic, combining two recent principles 
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of local search, Variable Neighborhood Search and Iterated Local Search. The main motivation 
behind the concept is the easy yet effective resolution of multi-objective optimization problems 
with an approach using only few parameters. 

After an initial introduction to the problem domain of flow shop scheduling under multiple ob- 
jectives, the introduced PILS algorithm has been applied to a set of scheduling benchmark instances 
taken from literature. We have been able to obtain encouraging results, despite the simplicity of 
the algorithmic approach. A comparison of the approximations of the Pareto sets has been given 
with a multi-operator local search approach, and, as a conclusion, PILS was able to lead to consis- 
tently better results. We had however to observe, that with increasing problem size, the number of 
iterations needed to identify an only locally optimal solutions grows exponentially. 
Nevertheless, the presented approach seems to be a promising tool for the effective resolution of 
multi-objective optimization problems. After first tests on problems from the domain of scheduling, 
the resolution behavior on problems from other areas might be an interesting direction for further 
research. 
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